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1 Introduction 

A class of non-equilibrium pattern formation problems appear when an interface 
moves with a velocity proportional to the gradient of some field u (which could 
correspond to temperature or impurity concentration in solidification problems, 
pressure or another appropriate potential in viscous fingering, etc.), which itself 
obeys a bulk equation (diffusion and Laplace equations respectively in these exam- 
ples) and a Dirichlet boundary condition on the moving interface. That constitutes 
a free boundary problem. Traditional methods to numerically treat these problems 
imply the explicit tracking of a sharp interface, whose dynamics is coupled to the 
bulk dynamics. This can be done either by treating both dynamics together with the 
corresponding moving boundary conditions or by projecting the whole dynamics 
into a single integrodifferential equation for the interface, explicitly nonlocal and 
highly nonlineal. 

The name phase— field model (PFM) or diffuse interface model denotes an al- 
ternative approach to the study of such problems. It can be seen as a mathematical 
tool that converts a moving boundary problem into a set of partial differential equa- 
tions, which allows for an easier numerical treatment. In these models an additional 
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scalar order parameter or phase field <p is introduced, which is continuous in space 
but takes distinct constant values in each phase. The physical interface is then lo- 
cated in the region where <fi changes its value, a transition layer of finite thickness 
W. The time evolution equation for <fi is then coupled with the u field in order to 
take into account the boundary conditions at the interface. When the equations are 
integrated the system is treated as a whole and no distinction is made between the 
interface and the bulk. The link with the original free boundary problem is ensured 
by requiring that it is recovered in the so-called sharp-interface limit W — > 0. 

In practice, the first phase-field models were derived in the context of solidifi- 
cation as Ginzburg-Landau time-dependent equations, dynamically minimising a 
free energy functional, in the spirit of model C of critical dynamics in the classifica- 
tion of Ref. [ 1 1. Most models for solidification-related phenomena continue to be 
derived from some thermodynamic potential. The sharp-interface limit then serves 
to relate the model parameters to those of the free boundary problem. Nevertheless, 
in other phenomena (see below) the existence of an appropiate thermodynamic po- 
tential might not be obvious, and the PFM equations have been constructed by 
directly "guessing" (see Sec. 13.21) what might reproduce the right free boundary 
problem in the sharp-interface limit, which is then completely crucial. From a 
computational point of view, such models are equally suited to simulate the origi- 
nal free boundary problem as long as they do reproduce it as W — > 1 . 

The earliest formulations of a PFM from model C were done by Fix J2J and 
Langer [3 1. A similar PFM was introduced by Collins and Levine [4 1, who applied 
it to the study of solidification from an undercooled melt with a diffuse interface, 
showing the existence of a solvability condition for the steady-state velocity of a 
flat interface for a continuous range of undercoolings. Moreover, they provided 
the first link between the PFM and the sharp-interface model with a kinetic term. 
Analytical properties of the equations in Langer's PFM have been studied by Cagi- 
nalp et al. ED ITU!, Gurtin and Fife and Gill OH- In Refs. I5IITT1. Caginalp 
et al. introduced anisotropic effects into a PFM. Caginalp also showed [71 El 
that various forms of the classical Stefan-type problem may be recovered as lim- 
iting cases of the PFM equations when W — > 0. Some critical situations where 
sharp-interface and PFM differ were presented in Ref. 1131 . Penrose and Fife fT4l 
provided a framework from which the PFM equations for the case of growth under 
non-isothermal conditions can be derived in a thermodynamically consistent man- 
ner from a single entropy functional. The earlier PFM, being obtained from a free 
energy functional, were only consistent in the isothermal case. 

Fix |2 introduced some numerical methods for the treatment of PFM equa- 

1 Indeed, some PFM for solidification cannot be completely derived from a single functional (non- 
variational models) either, but their discretization might converge faster (see Sec. l2.3> . 
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tions. After it, a large amount of numerical computations have been carried out. 
Early numerical computations in one dimension were done by various authors 
fT5l l81i. Kobayashi employed an anisotropic PFM for a pure substance in two [ 16 1 
and three [17| dimensions and simulated the evolution of dendritic-like structures 
growing into an undercooled melt. The results obtained showed the computational 
power of this new approach. 

Wheeler et al. developed and tested thermodynamically-consistent PFM for 
isothermal phase transitions in binary alloys ffSl l"i"9l and for the non-isothermal 
case in pure substances l20ll2T1l . Both models were used to analyze the linear sta- 
bility of a planar solidification front in Ref. J22), where a good agreement with the 
free boundary problem result was found when W was taken small enough. Sim- 
ilar results had been obtained by Kupferman et al. 1231 with a slightly different 
PFM. By combining aspects of the models in Ref. [ 18 1 and in Ref. EOl . Warren 
and Boettinger [ 24 1 constructed a PFM for the study of a non-isothermal binary 
alloy. The behaviour of the model derived in Ref. l20l was evaluated for varying 
W and spatial and temporal resolution [25]. This model was also used to calcu- 
late the operating states (tip velocity and radius of curvature) of dendrites growing 
from undercooled melts and to study the influence of anisotropic surface tension 
and kinetic coefficient [26]. Results obtained from this model were compared with 
solidification experiments in Ref. l27l . and this PFM was also used to study den- 
dritic growth in a channel l28ll . Moreover, it was applied to the study of mesophase 
growth in liquid crystals [ 29] |30l |^ |32) and it also permitted the study of growth 
with anisotropic heat diffusion I32ll33ll34l . 

Kupferman et al. 1351 used a PFM to study dendritic growth in the large under- 
cooling limit while a similar model was used by Mozos and Guo [36] in the limit 
of small undercooling. More recently, some properties of the solidification front in 
a supercooled liquid were derived using a PFM ll37ll . 

Fluctuations have been long introduced in the PFM in an ad hoc way, in order 
to destabilize morphologically unstable interfaces and obtain growing patterns with 
realistic characteristics, for example dendrites featuring sidebranching lIToTl . Noise 
in PFM can also account for external sources of fluctuations [32 38|. In a more 
rigorous approach, thermodynamical internal noise can be introduced in the PFM 
to account for equilibrium fluctuations l39ll40l . This approach allows the quanti- 
tative study of phenomena induced by internal fluctuations, like sidebranching in 
solidification of pure substances [39|. 

Since PFM results in principle depend on the ratio W/l c of the interface thick- 
ness to the smallest characteristic length scale of the original free boundary prob- 
lem l c as found in [25 1, simulations had to use very small values of it if agreement 
with the free boundary problem (where W = 0) within a few percent was desired. 
This became computationally prohibitive for experimentally relevant physical pa- 
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rameters where the system size is in turn much larger than l c , since the spatial dis- 
cretization must resolve W. Although physical interfaces are diffuse at the atomic 
scale, their thickness W cannot attain realistical small values even in simulations of 
mesoscopic structures, so that W j l c is necessarily exaggerated, and its effects can- 
not be disregarded. This large range of length scales to deal with has been the main 
drawback of the phase-field method. There have been two complementary ways 
of circumventing it: On the one hand, the explicit corrections to the original free 
boundary problem have been computed in a few cases, either as higher order terms 
in the sharp interface limit [^02] |43] |44l |45j or in the context of the so called thin 
interface limit l4^1l47ll48ll49ll50ll5Tl . Some of these corrections can be cancelled 
out I50ll43ll44ll45l . A complete cancellation has so far only been achieved for the 
solidification of a pure substance [46JE) or binary alloy into a single solid phase 
E8ll49l. It implies that the simulation outcome becomes roughly independent of 
W j l c at some finite value of it, and directly coincides with the original free bound- 
ary problem within a few percent. Quantitative simulations can then be run at that 
value, with no need to further decrease W/l c , which results in an enormous effi- 
ciency gain. Moreover, the possibility of such complete cancellation has the appeal 
of dealing with qualitatively distinct features such as zero kinetic term [ 46 1 or zero 
transport coefficient in one phase (one sided model) [48]. 

On the other hand, large system sizes can be treated exploiting the fact that 
the variation of the physical field u decreases with increasing distance to the inter- 
face. The resolution can hence be accordingly decreased, by either using adaptive 
meshes 1521 . or random walkers (whose amount naturally decreases) to simulate 
diffusion 1531 . This has enabled quantitative comparison between theory and ex- 
periments of dendritic solidification of a pure substance in three dimensions [54|, 
showing the power of the method. To fine tune the different physical anisotropics 
entering the problem, it is also necessary to control and eventually diminish the 
effects of the underlying grid anisotropy I35l l47ll55l . 

In an effort to further approach terrestrial solidification experiments in mas- 
sive samples, forced-fluid flow in two [56] and three [57] dimensions and natu- 
ral, buoyancy-driven convection [58| have only recently been incorporated in the 
phase-field formalism, usually requiring the adaptive mesh techniques mentioned 
above to solve the Navier-Stokes equations in the melt. PFM have also be formu- 
lated for two-phase flows [59|. 

Although the first phase-field models were intended for the solidification of a 
pure substance or binary alloy into a single solid phase, industrial alloys usually 
involve more components and / or solid phases. The first models for solidification 
into more than one solid phase focussed on eutectic growth, where two solids of 
different compositions can grow from a same melt, giving rise to composite struc- 
tures which are the most common solidification patterns after dendrites. Such first 
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models used the standard solid-liquid phase-field and the concentration field [ 60 1 
or a second, solid-solid phase-field [61 ] to distinguish between phases. Also using 
two different phase fields, peritectic growth has more recently been studied [62) . A 
formulation to treat an arbitrary number of phases using as many phase fields (in- 
terpreted as volume fractions) as phases was introduced by Steinbach et al. [63]. 
Later on, this formulation was used in combination with a free energy with cusp- 
like minima, which recovered for the first time pure two-phase interfaces as in the 
free boundary problem [64], removing the main difficulty to extend the progress in 
solidification into a single solid phase to model multi-phase growth quantitatively 
with a larger W/l c ratio. Since such pure two-phase interfaces have recently been 
also achieved for the three-phase case but using a smooth, non-singular free en- 
ergy [65], this possibility has become more real, and a refinement of this model 
has already allowed the first quantitative simulations in some situations ll50l . 

The phase-field method is also quickly expanding in different directions: Be- 
yond applications in solidification (for a review on some of the topics discussed so 
far, see Ref. 1661 . with a mention in solute trapping llT9ll67l ). PFM have been de- 
rived for other phase transitions and transformations as those between mesophases 
in liquid crystals mentioned above |29l|30||21, solid-state structural phase trans- 
formations, and electrodeposition. Other PFM with an impact on materials science 
include nucleation [68 1, grain growth [69| and coarsening, domain evolution in thin 
films, spiral growth [ 70 1, and models incorporating elastic effects, such as those for 
stress-induced instabilities [71], the dynamics of dislocations [72 1 and the propa- 
gation of cracks [73]. (For a short review on some of those, see Ref. ]74|). 

Further outstretching their reach, some PFM have recently been derived for the 
dynamics of fluid interfaces. Because such cases are examples of interfaces be- 
tween media which are not necessarily different thermodynamic phases of a same 
substance, the existence of a free energy functional to minimize is not obvious. 
The first model for viscous fingering including arbitrary fluid viscosities, for in- 
stance, was derived purely on the basis of reproducing the desired free boundary 
problem in the sharp-interface limit, without the guide of any free energy func- 
tional ll43"l l44l. In contrast, later formulations of the same problem were derived 
from such a functional [75 1. This nicely illustrates two crucial points: On the one 
hand, given a free boundary problem, a PFM for it is not unique, which allows for 
the necessary freedom to construct models cancelling some of the finite interface 
thickness corrections |4l] |42l |43] @4l |45] |46l |47) to that free boundary problem, 
yielding more accurate and computationally efficient models, as discussed above. 
On the other hand, a PFM can derive or not from a free energy functional. We 
emphasize this latter point, because we believe that removing this restriction might 
help derive models for other problems. An even more recent example of that is 
a PFM for vesicle dynamics [45|, based on the viscous-fingering model of Refs. 
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I43ll44l and therefore not derived from a free energy functional. 

Finally, let us mention that the key step of the PFM, the calculation which 
relates it to a specific free boundary problem, can sometimes also be applied to 
phenomena with a physical diffuse interface in order to obtain an effective sharp- 
interface description of them for theoretical purposes. One example of those are 
the steady states of thermal plumes, where finite interface thickness corrections to 
the resulting free boundary problem explained the selection mechanism for their 
width [42|. Not all diffuse-interface problems can be mapped to an effective free 
boundary problem. It has been shown that the so-called pulled fronts cannot [76], 
and other singular fronts require a generalisation of the method employed [771, as 
those arising in reaction-diffusion systems with a non-linear diffusivity, e.g. in 
physical models of bacterial colony growth 1 771 . 

The main features of PFM that make them adequate to describe interfacial dy- 
namics can be summarized as follows: 

• As sets of coupled partial differential equations, they are simpler to inte- 
grate numerically than the integrodifferential type of equation usual in sharp- 
interface descriptions. (In particular, they do not suffer from long-time nu- 
merical instabilities as some of the latter). This also makes modifications 
of the model comparatively easier to incorporate, e.g. to include anisotropy 
EUEEZlEa, external fluctuations HUE! or fluid flow EH EU EU . 

• PFM are especially suited for the introduction of internal fluctuations, since 
these are usually known from first principles [39]. 

• More specifically, PFM have the well-recognised advantadge that the loca- 
tion of the interface does not have to be explicitly determined as part of the 
solution, but it is obtained from the solution of the equations for the phase 
field(s). Not having to track the interfaces, topological changes of the in- 
terface (self-intersections, coalescence or pinching of droplets, etc.) and 
extensions to three dimensions are more feasible within the PFM formalism. 
(Sharp-interface formulations are normally restricted to connected interfaces 
and to two dimensions, either by the formulation itself or because they be- 
come typically impractical in those other cases). 

Recently, the concept of PFM has been used in a broader sense to include 
any model which contains continuous fields that are introduced to describe the 
phases present with diffuse interfaces I5T1 l78l l79l . In this context models for 
order-disorder transitions, spinodal decomposition and heterogeneous nucleation 
II5TII79I . described by only one equation for a scalar field, has been included at the 
same level than dendritic growth and solidification of eutectic alloys. In Ref. ll78l 
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a continuous one-sided model has been proposed to model Hele-Shaw flows in the 
high viscosity contrast regime. In the same spirit, continuous models have been 
introduced in the field of roughening due to their nonlocal properties [801 I8T1 . 

In this chapter we review some results concerning the use of PFM in the context 
of pattern formation in morphologically unstable interfaces. More specifically, we 
address two problems which have played a prototypical role in this context, namely 
solidification and viscous fingering. Whereas PFM for solidification are usually 
derived from a free energy functional, this is not necessarily the case for viscous 
fingering, as mentioned above. This will enable us to illustrate the two situations: 
In Sec. |2]we derive a PFM for solidification from a free energy functional; in Sec. 
|3] we discuss the appropriate PFM for viscous fingering which directly gives the 
sharp interface model. In both cases the asymptotic calculation that relates the 
PFM to the sharp interface counterpart is a central step to validate the model. In 
the next sections we will show as examples results of applying PFM to different 
situations: fluctuations in solidification, growth of liquid crystal mesophases which 
are characterized by strong anisotropies, and viscous fingering in Hele-Shaw cells. 



2 A first example: solidification of a pure substance 
2.1 The Sharp Interface Model 

The standard sharp-interface macroscopic model for free solidification of a pure 
substance [ 82 ] considers the solid-liquid interface as a microscopically thin mov- 
ing surface. This model relies on the heat diffusion equation together with two 
boundary conditions at the interface, which is assumed to be sharp. The diffusion 
equation for the temperature field T is given by 

BT 
at 

where t is time and D is the diffusion coefficient (D = k/c p , being k the heat 
conductivity and c p the specific heat per unit volume). Tilded variables denote 
magnitudes in physical units. We are considering the symmetrical model of solidi- 
fication, that is, we assume that D and c p are equal in both phases. 

One boundary condition is obtained from the conservation of the released latent 
heat at the moving interface: 

Lv n = Dc p [(V n T) s - (V„T) L ], (2) 

where L is the latent heat per unit volume, v n is the normal velocity of the inter- 
face, V n is the normal derivative at the interface and S and L refer to solid and 
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liquid respectively. The left-hand-side corresponds to the rate at which heat is 
produced at the interface per unit area while the right-hand-side is the total energy 
flux away from the interface. The sharp interface model is completed with a ther- 
modynamic boundary condition which assumes local equilibrium at the interface 
and incorporates the non-equilibrium effect of the kinetic attachment of atoms at 
the moving interface: 

T M cr _ v n 

J-int = J-M y W , (->) 

L \i 

where Ti nt is the temperature at the interface, a is the surface tension, fj, is the 
kinetic coefficient and k is the local curvature of the interface, being positive when 
the solid bulges into the liquid. This relation can be thermodynamically derived 
by considering that at local equilibrium the solid and liquid Helmholtz free en- 
ergies are equal. The ratio a/L has dimensions of a length which sets the scale 
of the pattern. If a = the so called Stefan problem is recovered and no real- 
istic pattern formation can be described. Introducing the reduced diffusion field 
u = (T — Tm) / {L / c p ) and scaling Eqs. (171131 with the diffusion length I = D/v c 
and the diffusion time 7 = l 2 /D (v c is some characteristic front velocity, in 
most experiments in the range of 1/xm/s), we obtain a dimensionless set of sharp- 
interface equations 

Vn=[(V n u) S -{V n u) L ], (5) 

I -ill 

where do = Tm<Jc p /L 2 is the capillary length of the substance and (3 = c v j (fJ,L). 
2.2 Derivation from a free energy functional 

PFM consist in a set of equations describing the dynamics of a continuous or- 
der parameter <p(r, t) which takes different constant values in the solid and liq- 
uid phases. The solidification front is then determined by the implicit condition 
4>(r, t) = const. , and the transition between solid and liquid phases takes place in 
a diffuse region of thickness W. The dynamics of the phase-field 0(r, t) is cou- 
pled with the evolution of the diffusion field, and can be derived from a free-energy 
functional such as 



F[e, 0] = J dv [K\ V0| 2 + h f(4>) + e u< 



(7) 



where u = e/eo + h(4>) /2, being e the local enthalpy per unit volume, eo = , 
and h a certain polynomial function of <p. The gradient term in the functional 
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takes into account the free energy cost of interfacial deformations and /(</>) is a 
double- well potential describing the bulk free-energy density. The minima of /(</>) 
determine the constant values of <\> at the bulk phases. The constants K and ho have 
dimensions of energy per unit length and energy per unit volume respectively. 

The basic phase-field equations can be then obtained by applying the well 
known variational relations 

d~ t e = De V 2 (jpj - V-q e (r,t). (9) 

These are the equations of model C, applicable to the dynamics of two fields, 
one conserved and the other nonconserved [ 1 ]. The terms fj(r , t) and q e (r, t) are 
noises accounting for thermal fluctuations. These noise terms can be ignored if 
fluctuations are not relevant in the study at hand. 

After some calculations and scaling space and times with 1, 7, we get the set of 
equations 



a£ 2 d t <P = e 2 V 2 4>- f'(4>) - ^-ti(<f>)u + V (r,t) (10) 

ho 

1 



d t u = V 2 u + -d t h- V-q u (r,i) (11) 

where e = ^f, W = ^jKjho and a = ■ The statistical properties of the 
noises can be determined by considering the fluctuation-dissipation theorem in Eqs. 

( V (r,t) V (r' ,t')) = ^0Ls(r-r')5(t-t') (12) 
<g*(r,t)g£(r' ) 0> = ^^^6 tJ 5(r - r')5(t - t'). (13) 

2.3 The thin-interface limit 

In this subsection we will perform the W — > limit of a PFM for solidification. 
Our aim is first to obtain the relationship of the PFM parameters with those of the 
sharp interface model, by using the thin interface asymptotics. Second we want to 
illustrate how can one construct a PFM with the only criterium of reproducing a 
desired sharp interface model. In particular we will use a variant of Eqs. (110111b 
derived above, not necessarily obtainable from any energy functional. 
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Consider the general model equations 



ae 2 d t 4> = e 2 V 2 </> - /'(<£) - Xg'(4>)u (14) 
d t u = V 2 u + ^d t h, (15) 

where A = |^ and no particular choice for the function h has been taken. For 
the particular choice of h(<p) = %g((j>) with b = g(l) — g(— 1), Eqs. (fl4l 1731 
recover the variational structure described in the last section. Other choices for the 
function h improve the computational efficiency of the model but the variational 
properties of the equations are lost. The non- variational formulation of the problem 
has some thermodynamical drawbacks: For instance, the fluctuation-dissipation 
theorem loses its validity and the noise properties should be calculated by means 
of other procedures [83]. In any case, the only restriction which must be imposed 
to the function h is that it satisfies h(±l) = =Fl. Usual choices for the free-energy 
density potential / and for the coupling function g are 

/w - (i« 
m (17) 

We present the thin-interface limit without considering the noise terms (the exten- 
sion to fluctuating PFM will be presented elsewhere [ 831). 

As we mentioned in the introduction, the phase-field Eqs. (11011 lb reproduce 
the sharp-interface problem when the interface thickness W is small enough com- 
pared with the smallest physical length scale, which in solidification is the capil- 
lary length do. However, Karma and Rappel introduced the so-called thin-interface 
limit, in which they assumed e = j € 1, but not necessarily W < do [47]. This 
corresponds to take the limit e — » while keeping A, which we will see to con- 
trol the ratio W/do, constant. Of course, it introduces corrections at first order in 
W/do coming from the term —\g'{<f))u. The point is that, for h{4>), g{<j>) and /(</>) 
odd, these corrections can be tracked to change only the value of (3 being simu- 
lated. Once one has identified the new (3, the model parameters can be adjusted to 
reproduce the desired j3 while keeping W even larger than do in practice 2 . This 
is completely equivalent to perform a more traditional asymptotic expansion using 
W/do as a small parameter but keeping corrections up to first order in W/do as 
proposed in the introduction. Indeed, corrections to the mass conservation Eq. ©, 
not coming from the —Xg'(<p)u, are ignored in the W/l — > limit. They turn out to 

2 As a side benefit, /3 can now be made to vanish, because these corrections can have the opposite 
sign that the W/do — > contribution. 
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vanish for h(<p) odd in the two-sided case considered here where D is the same in 
both phases, but not in the one-sided case where the solid diffusivity is negligible 
compared to that of the liquid I4fll48ll49l . Here, we follow the thin-interface limit 
in Ref. 1371 . as the most compact way of obtaining the correct result. For a higher 
order W/do expansion, see the appendix of Ref. [47], [41 1 (for unequal diffusivi- 
ties) or [49 1 (for the strict one-sided case) for more detailed presentations. See also 
Ref. I5T1 . An expansion without any a priori assumption on which physical scale 
to compare the phase-field small parameters with can also be found in Ref. [43 1 in 
the context of viscous fingering. 

In order to perform the asymptotic expansion, we divide our system in two 
different regions: an outer region at the bulk phases far away from the interface 
and an inner region located around the interface at a distance W. We will use capital 
letters to refer to all the fields in the inner region. In the outer region we expand 
both fields for e <C 1 

u = n + eui + O(e 2 ) (18) 
<f> = <J)o + efa + 0(e 2 ). (19) 

Introducing Eq. (fl9l into Eq. H4i we obtain, at order zero in e 

-f'(<h) ~ A</(0oH = 0. (20) 

It is easy to see that <po = ±1 is a solution of the last equation for our particular 
choice of the functions /(</>) and g{4>). As h is constant in the bulk, the outer 
solution of the diffusion field at zero order is given by the diffusion equation 

d t u = V 2 n . (21) 

Higher orders in e are described by 

<k = 0, (22) 
d tUi = V 2 Ui , i > 0. (23) 

In order to work in the inner region, we consider Eqs. (fl4l and dT3l and write them 
in a three dimensional orthogonal curvilinear coordinate system (r, s, to) based in 
the surface defined by </>(r, t) = 0. The normal coordinate of this coordinate system 
r is then scaled using the small parameter e so that p = |. Writing the equations 
in the frame of the moving front and keeping terms until first order in e we have 

aved p $ = d 2 p <S> - f + £Kd p $ - Xg'U, (24) 
-£vd p U = d 2 p U + eKd p U - y3 p /i($), (25) 
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where k = k s + k w being k s and k w the two principal curvatures of the level-set 
<fi = 0. Inserting the inner expansion 

$ = $ + e $ 1 + o( e 2 ) (26) 
[/ = U Q + eC/i + 0(e 2 ) (27) 

into Eqs. (l2"4l and d23l we get, at zero order 

a p 2 <I>o - /'($ ) - A 5 '(<Wo = (28) 
d 2 [/ = 0. (29) 

A valid solution of Eqs. d28ll29T ) is given by a kink for the phase-field 

$o(p) = -tanh(-^), (30) 

and a vanishing u at zero order, i.e. Uq = 0. At first order we get 

J2$! = -(va + K)d p $ + \J($ )Ui (31) 
d 2 ^ - |^($ ) = 0, (32) 

where 17 is the linear self-adjoint operator given by 

n = d 2 p - /"(*„). 03) 

Integrating Eq. d32l over p we obtain 

d p U 1 = A + ^h($ ), (34) 

where A is an integration constant independent of p. Integrating Eq. (13*41 over p 
we get an expression for the diffusion field at first order 

U 1 = B + Ap + ^- f dp'h($ ), (35) 
^ Jo 

where B is another integration constant. To determine the value of the constant B 
we impose the solvability condition for the existence of a solution for $x, which 
can be written as 

/oo 
dp (d p <3? ) = (36) 
-oo 

and allow us to find an expression for B 

Jva + K)h vh 
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where the integrals I\ , I2 and ^3 depend on the particular choice of the functions g 
an h and are given by 

/oo 2\/2 
dp(d p $ ) 2 = (38) 
-00 <j 

/oo 
^(a p * y(*o), (39) 
-00 

dp(5 p * y(*o), / <W$ )- (40) 

-00 Jo 

In order to determine the value of the integration constant A, we impose the match- 
ing of the inner and outer diffusion gradients at the interface. This procedure will 
bring to an expression for the conservation of the diffusion field across the inter- 
face. The matching conditions for the derivatives of u can be written in the form 

lim d r u(r,t) = lim -d p U(p,s,w,t). (41) 

r-»0± p-^±oo £ 

We insert then the e expansions into Eq. d4TT > and equal the corresponding orders. 
For the matching of the derivatives at zero order we get 

d r u \ ± = A±^, (42) 

where we have used Eq. d34b and that h(&o(p — > ±00)) = =f1, which is the only 
restriction we have required for the function h. This last equation can be used to 
find a value for the constant A 

A = ± V - + d r u a \ ± . (43) 

This last equation can be used to get an expression for the conservation of the 
diffusive field across the interface 

v = <9 r uo|~ — d r uo\ + , (44) 

which is the form of the heat flux conservation across the interface appearing in the 
sharp-interface problem Eq. 

The last thing we have to do is to recover the Gibbs-Thomson equation in the 
thin-interface limit of the model. For this we impose the matching conditions for 
the diffusion field 

lim u(r,t) = lim U(p,s,w,t), (45) 

r—>0 ± p^ioo 



(46) 



and using Eqs. d351 1371 and (1431 we obtain 

"aii I3 I4 
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where the integral I4 is given by 



±00 



dp [fc(*o)±l]- 



(47) 



Equation (l4*6l l can be directly compared with the Gibbs-Thomson equation in the 
sharp-interface description (Eq. |6jl. Identifying terms we can get an expression for 
the phase-field parameters in terms of the sharp-interface physical constants 



da 



ali 



2/2 



h 

2 



(48) 
(49) 



Note that in the last equations the integrals I 2 — I4 depend on the particular choice 
for the functions g and h. Another interesting feature of this result is that the 
expression for j3 allows to simulate situations with neglecting kinetic attachment 
by means of an appropriate selection of the phase-field parameters. 



3 A second example: The Saff man-Taylor problem. 

The so-called Saff man-Taylor or viscous fingering problem has played a central 
role in the context of interfacial pattern formation because of its relative simplicity 
both experimentally and in its theoretical formulation |84||85l|86l|87l|88l|l9l. It 
deals with the destabilization of the interface between two immiscible fluids when 
a less viscous fluid is injected to displace a more viscous one, and / or when a more 
dense fluid is placed on top of a less dense one in a two-dimensional geometry 
known as a Hele-Shaw cell. This consists of two parallel glass plates of dimensions 
L x x L y separated a distance b <C L x , L y . 

In contrast with the three-dimensional Rayleigh-Taylor instability, in such a 
cell and for some range of injection velocities Voo, the flow is overdamped (small 
Reynolds number limit). Therefore, the velocity u is proportional to the gradient of 
a potential defined as the deviation from hydrostatic pressure: u = — (b 2 /12fi) V (p— 
P9csy), where p is the pressure, \i the shear viscosity, p the density and g c g an ef- 
fective gravitational acceleration in the y direction for non-horizontal cells. 



3.1 The sharp-interface model 

Assuming incompressibility in the bulk, V • u = 0, continuity of normal velocities 
and a pressure drop proportional to the interfacial tension a on the interface, one 
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obtains the classical free boundary problem 



V 2 p = 0, (50) 



(Vpi + pig c sy = Vp2 + /025efry = -To v m ( 51 ) 

PI-P2 = (TK, (52) 



where the first equation holds in the bulk of each fluid, and the next are its boundary 
conditions on the interface between them; subscripts 1 and 2 label each fluid and 
the notation pi stands for the pressure on the interface coming from fluid i; r is a 
coordinate normal to the interface, v n , its normal velocity, and k, its curvature. 
The dynamics are controlled by the two dimensionless parameters 

B = ^T2im ^ 7 ^>°< c=^^>0. (53) 

nLUVooim - /j, 2 ) + g e ff{pi - P2)\ + 

B is a dimensionless surface tension, and can be understood as the ratio between 
the capillary (stabilizing) force and the driving (destabilizing) force (injection+ 
gravity), and c is the viscosity contrast, which is so far completely arbitrary: < 
c < 1. This corresponds to having set ourselves in the frame moving with the fluid 
at infinity (or, equivalently, with the mean interface) and taken L x as unit length 
and [/* = cVoo + g ^ 2 as unit velocity (see Ref. (90)). 

Analytical understanding of these highly nonlinear and nonlocal dynamics [Eqs. 
d50H52l l is basically restricted to high viscosity contrast c = 1 and small surface 
tension (9TJ |92j |93]| , so one relies mostly on numerical work I901l93ll94ll95ll96ll971 
l98l . A systematic weekly nonlinear analysis is also possible for the early stages of 
the evolution and arbitrary parameters (99). In the next subsection we introduce a 
PFM for Hele-Shaw flows with arbitrary viscosity contrast c. Although in the high 
contrast limit c = 1 the Hele-Shaw dynamics are quite analogous to the one-sided 
solidification problem (in the appropriate approximations [96]), the arbitrary vis- 
cosity contrast case has been shown to exhibit quite different dynamics, and has 
raised some interesting questions, particularly concerning the sensitivity of finger 
competition to viscosity contrast I90l l94ll95lll00 l and the long time asymptotics 
of the low viscosity contrast limit c = I95lll01l . 



3.2 Discussion of the sharp interface limit 

As emphasized in the introduction, the classical Saff man-Taylor problem does not 
involve phase transformations, so that a derivation of a PFM for it from a free 
energy functional is not necessary. Although such a variational model is possible 
[75], here we discuss the first PFM for the Saff man-Taylor problem with arbitrary 
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viscosity contrast of Ref. l43l . without the help of such a functional, as done in that 
reference, in order to illustrate how such a derivation (and the subsequent sharp- 
interface limit) work for a particular case. The idea is to construct two coupled 
partial differential equations, one for some scalar field describing the flow and 
one for a phase field, such that we recover Eqs. ( l5T)H52t in the sharp-interface 
limit. How this can be achieved depends of course of the particular free boundary 
problem. 

An important difference between the arbitrary c Saff man-Taylor problem Eqs. 
d50ll52l and the solidification one, as stated by Eqs. ([BGll, lies in their boundary 
conditions: in the two-sided solidification problem the physical field u is continu- 
ous and its normal gradients, discontinuous, through the interface; in the Saffman- 
Taylor problem it is the other way round, a velocity potential is discontinuous and 
its normal gradients are continuous (because the velocities are so). We therefore 
begin by rewriting the free boundary problem in terms of the harmonic conjugate of 
the velocity potential, the stream function ifi. This exchanges normal and tangential 
derivatives (x • u = d y xp, y ■ u = —d x xp) so that we obtain, in the adimensionalisa- 
tion mentioned in the last section, 

VV = 0, (54) 

Vv(0+) -Vv(O-) = -7-c[</y(0+)+Vv(0-)], (55) 
^ s (0+) = ^(0-) = -v n , (56) 

where s is arclength along the interface, subscripts stand for partial derivatives 
except for v n , and 

= Bk s + y ■ s, (57) 

with k(s) the interface curvature. 

As we can see, the boundary conditions continue to be substantially different 
from those of two-sided solidification, since the stream function rp is now con- 
tinuous, but the value it takes on the interface and the magnitude of its normal 
derivative jump are not the same than in solidification. Nevertheless, because ift 
is continuous through the interface, a single field (and not one field in each fluid) 
can be used for the whole bulk, which is obviously advantageous when it comes to 
build a PFM, since we recall that PFM treat all the system as bulk. 

Note that Eqs. (I54l55t can be written together as 

V 2 ^ = -w, w = { 7 (s) +c[t/v(0 + ) +^ r (<r)]}$(r) (58) 

where 5(r) is the Dirac delta distribution and w = z ■ (V x u) is the fluid vorticity, 
which is confined to the interface. This suggests to construct our PFM equation 
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for ip as a regularisation of this distribution equality, coupled to the remaining 
boundary condition, the continuity of normal velocities, by advecting the phase 
field 9 (and therefore the interface) precisely with the physical fluid velocity on the 
interface. 

More specifically, the equations proposed in Ref. [43] read 

dib -> - 11 

€ ~it = V V + cV -( 0V VO + -^7(#)(l-O, (59) 

89 

e 2 — = f(9) + e 2 V 2 9 + e 2 K{9)\V9\+e 2 z-(V^xV9) (60) 

where f(0) = 9(1 - 6 2 ), and ^f - = s{9) ■ (BVk(9) + y), «(0) = -V • f(9), 
with f(9) = and s(9) = r(9) x z. j(9), k(6) are functionals which generalize 
the magnitudes defined above for the interface to any level-set of the phase-field, 
and e is a small parameter of the model; in particular it can represent the interface 
thickness W (see below). This model is inspired in the vortex-sheet formulation of 
the problem [90]. Similar ideas have previously been applied to describe physically 
diffuse interfaces in the context of steady state selection in thermal plumes [42]. 

If we leave the two last terms aside, Eq. d60l) is the time-dependent Ginzburg- 
Landau equation for a non-conserved order parameter or model A (without noise) 
in the classification of Ref. [ 1 1. The field in this model is known to relax towards 
a kink solution in a short time scale, and then to evolve to minimize the length 
of the effective interface according to the Allen-Cahn law (i.e. with normal ve- 
locity proportional to the local curvature). The factor multiplying the laplacian 
has been chosen to be e 2 for the kink width to be 0(e), so that e can be consid- 
ered the interface thickness. On the other hand, the e 2 factor in the time derivative 
ensures that the relaxation towards the kink is much faster than the evolution of 
the interface. Notice that model A describes the relaxational dynamics of a non- 
conserved order parameter, whereas our problem is actually non-relaxational and 
strictly conserved (mass conservation and immiscibility). The other two terms in 
the phase-field equation will correct this apparent contradiction. In order to cancel 
out the local Allen-Cahn dynamics of the interface which is built in model A, we 
add the term e 2 k(9)\V9\. Such a term cancels out Allen-Cahn law by giving rise, 
to leading order, to an identical contribution but with opposite sign. 

With these elements so far, our phase-field relaxes to a kink profile located 
along an arbitrary interface which remains almost completely stationary, regardless 
of its shape, provided that the interface thickness remains smaller than the local 
radius of curvature (i.e., that a sharp-interface description makes sense). This is 
because the dynamical effect of surface tension associated to the Ginzburg-Landau 
free energy of model A without noise has been removed (up to first order) and 
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the interface has not yet been coupled to the fluid flow, represented by the stream 
function. This coupling is achieved by adding the last term in Eq. d5Ul . which 
stands for — e 2 u • V# and thus sets the phase-field — and therefore the interface — 
in the frame moving with the fluid velocity u. This term restores the fully nonlocal 
dynamics of the Hele-Shaw model. In particular it yields the continuity of normal 
velocities Eq. d56l and reintroduces surface tension, which is contained in the 
dynamical equation for the stream function through 7(6*). 

At first sight, the idea presented in Ref. [43 1 and reviewed here of cancelling 
out some undesired effect (here, the Allen-Cahn law) from the PFM equations by 
adding an extra term designed to produce the same effect but of opposite sign close 
to the sharp-interface limit [here, the term e 2 K(9)\V9\] might seem somewhat arti- 
ficial or "unphysical" compared to PFM of solidification where the surface tension 
arises more "naturally" from the free energy cost for the gradients of the phase field 
without the need to cancel it out and reintroduce it in a different way. Neverthe- 
less, it has been successful in later developments, both in the sharp-interface limit 
and at first order in the interface thickness. The present formulation has recently 
been adapted to treat vesicles [45 1. Seen as an interface, the vesicle membrane not 
only evolves to conserve mass inside and outside the vesicule, but also keeps its 
own length roughly constant. Consequently, another term was added to cancel out 
interface length changes, yielding a suitable model to treat such membranes [45|. 
Similarly, an extra term in the standard PFM of solidification cancels out a spurious 
solute trapping correction to the the one-sided model at first order in the interface 
thickness J48]|49|. Combined with an appropriate choice of the model functions, 
all first-order corrections can then be cancelled out ll48ll49l. 

As for Eq. d59t . its right hand side is intended to reproduce Eq. d58l . and 
therefore also Eqs. d54l and (I55t . If the phase-field 9 has a kink shape, 1 — 8 2 is 
a peaked function which, when divided by e, gives rise to the delta distribution for 
the vorticity. However, this only accounts for the 7 in the weight of the delta. The 
part proportional to the viscosity contrast c must be introduced separately as the 
cV • (9Vip) term because of the non-local character of ip r (0 + ) + ip r (0~). Finally, 
the time derivative is multiplied by e to recover the laplacian (and not diffusive) 
behavior of Hele-Shaw flow in the sharp-interface limit. 

In the PFM the interface width and the convergence to the sharp interface limit 
is controlled by the small but finite value of the parameter e. Note that the role of 
the sharp-interface limit in this PFM is not the same than in solidification, in the 
sense that the parameters of the model do not need to be identified, but are just 
"built in" in the model. The sharp-interface limit hence merely serves as a safety 
check. Nevertheless, one can still wonder which value of e is needed to accurately 
reproduce the actual Hele-Shaw dynamics for given values of the physical param- 
eters B and c. This question can be qualitatively answered by noting the distinct 
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roles played by e in the phase-field equations, Eqs. (I59l60t : 

The e factors appearing in e 2 V 2 6, e 2 K(8)\V6\ and - 2 ) all stand 

for the interface thickness, and this is required to be small compared to the longi- 
tudinal length scale (A;^ 1 of the interface: e|A;| << 1. 

In contrast, the e in e|^ has nothing to do with the interface thickness (and we 
will therefore denote it by e from now on), but its aim is to ensure that the stream 
function is laplacian and not diffusive in the e — > limit, which commutes with the 
e — > one: e sets the time scale of the diffusion of the stream function through a 
given characteristic length of wavenumber k, ^A^ , which must be much smaller 
than the characteristic growth rate of the interface , so that the stream function 
is slaved to the interface: 5& << 1 ± c. 

The e 2 in e 2 ^f represents the relaxation time of the phase field towards the 
steady kink solution (see Eq. t60b . which must be kept well below the interface 
growth time l^l -1 for the phase-field to remain close to the kink profile during 
the interface evolution: e 2 \cu\ << 1. This factor must be the same that the one in 
e 2 z ■ (Vtp x W) in order to recover the macroscopic equation Eq. d56l . In fact 
there are at least two distinct powers of e for this relaxation time (e and e 2 ) for 
which the right sharp-interface limit is achieved, and the corrections to it are also 
the same. 

To sum up, there are at least two independent small parameters (e and e) con- 
trolling the limit. When trying to approach macroscopic solutions by means of 
numerical integration of the phase-field equations, it is very convenient to vary 
them independently in order to save computing time, since both affect it [44]. 

A more quantitative answer to the question of the necessary values of e, e 
to achieve a given precision can be obtained by computing not only the sharp- 
interface limit, but also the corrections to it at first order in the interface thickness 
e considering e of 0(e). Thus, one obtains an effective free boundary problem 
(where the interface is mathematically sharp) keeping track of corrections partly 
due to the fact that the interface of the PFM is not sharp. More precisely, it keeps 
track of corrections up to first order in both e and e. This limit follows basically the 
same procedure described in Sec. 12.31 It is performed in Ref. 1431 . where a detailed 
account of how the successive orders in the asymptotic expansions are derived and 
matched can be found. Here, we just give the result for the sake of brevity: 
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= - Vn + eV2^Tc+^c 2 + 0(c 3 )}, (63) 

ZD 5 

where T = 7 + c[ip r (0 + ) + ip r (0~)] is the weight of the vorticity defined in Eq. 
d5£t evaluated up to 0(e) and 

9± (c) = l-^ + ( ± l + i^i)lni±^. (64, 

This contains two important results: On the one hand, the correct free boundary 
problem of Eqs. (15415 61) is recovered in the sharp-interface limit e,e — > 0, as 
expected. On the other hand, the corrections at 0(e, e) in Eqs. doTT) and d63b are 
obtained and go as e and e respectively, whereas Eq. d55b remains unaffected. Note 
that the correction in e appearing in Eq. d63l has nothing to do with an Allen- 
Cahn law. So the e 2 K(9)\V9\ term has cancelled that effect out even in the first 
order corrections. This is the partial improvement of the convergence to the sharp- 
interface limit achieved by this model and mentioned in the introduction. 

4 Applications 

In this Section we present results obtained in simulations of PFM in different pat- 
tern forming systems. We will start with two examples taken from solidification 
problems, namely directional solidification of alloys, and dendritic sidebranching 
in free solidification. In these problems we will devote special attention to the ef- 
fect of fluctuations. Next we will apply the same kind of models to the growth 
of mesophases in liquid crystal. This problem is mathematically similar to those 
of solidification, but with the particularity of the presence of strong anisotropics, 
both in the interface and in the bulk. We will show how the PFM is able to deal 
these anisotropics and to reproduce the morphological transitions observed in ex- 
periments. Finally, we present some results for viscous fingering obtained with the 
PFM discussed above. 

4.1 Directional solidification 

The first application will be the study of initial transients in directional solidifica- 
tion of an alloy in the context of the symmetric model ll82l . This example will be 
used to show the numerical convergence of the PFM to the sharp interface model. 
Once the convergence is reached, the PFM will be used to quantitative test of the- 
oretical predictions. This will be done in the context of the symmetric model of 
solidification. To that end we will present predictions for the transient recoil of the 
sharp interface model and for the transient dispersion relation by using an adiabatic 
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approximation. We will then perform phase field simulations of this initial regime, 
both for the evolution of single modes and for the complete system with thermal 
fluctuations. 

The selection of a dendritic pattern during the directional solidification of a di- 
lute binary alloy is a complex problem which depends on initial conditions! 102 ], 
and in particular on the first wavelengths that appear in the destabilization of the 
planar front induced by fluctuations [ 1031 11041 . Recent work has focussed in the 
importance of internal fluctuations in solidification patterns 11051 11061 . Quanti- 
tative agreement with experiments has only been obtained for the solidification of 
pure substances [ 107 ], whereas in solutal growth the origin itself of the fluctuations 
is still an open problem [38 1. The analysis of the evolution of single modes and 
of fronts with fluctuations, both from theory and from simulations, should be of 
relevance in the dendritic selection problem. 

4.1.1 Predictions of the sharp-interface model 

In a directional solidification experiment, a thermal gradient Gz is moved in the z 
direction along the sample at constant pulling velocity Vz. Provided the sample is 
thin, we take the system as 2D and describe the interface position in the moving 
frame of the gradient by z = £(p,t), where p = xx + yy. The solid phase is 
located in the region where z < £, and the liquid where z > £. We will consider 
the particular case of symmetric directional solidification, which assumes the same 
solute diffusivity D in both phases. We introduce diffusion length I = D/V and 
time r = D/V 2 to scale variables as r = f/l and t = t/r. We also introduce a 
diffusive field in each phase Ui(r,t) = C ^^//° (i=l solid, i=2 liquid), where Cj 
is the solute concentration, the concentration far away from the solid-liquid 
interface, and ACo = [C 2 — Ci]i n t the concentration jump across the interface. In 
the moving frame and in reduced variables, the fields Uj evolve according to the 
diffusion equation 

(|-|-V>, (r - t ) = 0. (65) 

The diffusion fields Uj must satisfy some moving boundary conditions at the 
interface position which can be written as 

[til - u 2 ]int = -1 (66) 
u 2 \ int = l-^ + fK (67) 

h-[V Ul -Vu 2 ] mt = -n z (l + i) (68) 

Eqs. d65W68t define the so called sharp-interface description of the symmetric 
directional solidification problem. Eq. d66l relates solute concentrations at both 
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sides of the interface. We have assumed the additional approximation of having 
a constant concentration jump across the interface. This particular assumption is 
equivalent to suppose that the mixture has parallel solid and liquid branches in the 
T{C) coexistence diagram, which is valid only for liquid crystals and alloys with 
a partition coefficient close to 1. Eq. do7l is the Gibbs-Thomson equation (local 
equilibrium at the interface). In this relation It = ACo/G is a thermal length 
imposed by the temperature gradient, and do is the capillary length. Eq. d68l 
describes the solute conservation across the interface, and h represents a normal 
unitary vector pointing to the liquid. 

Using Green's function techniques [ 108 1, it is possible to derive a closed in- 
tegral expression for Ui at each side of the interface. Introducing the notation 
p = (p, z,t), ps = (p, £, t), the corresponding integro-differential equation reads 

Jj dt'Jdp'(l + i(ff,t))G(p s ,p' s ) + i(l + erf{ * m -MM +t )) 

where G(p,p') is the Green function for Eq. d65l . and £ = stands for the steady 
position of an unperturbed planar interface. Note that Eq. do^t includes transients 
from the initial condition at to. The next step is to perform a linear stability anal- 
ysis of the problem to obtain a transient dispersion relation describing the time 
evolution of a sinusoidal modulation with wavevector kl. Within an adiabatic ap- 
proximation 11031 11041 we derive the growth ratio of the mode as 

u(kl, t) = ±£ _ (1 + £) 2 + q 2 (l + - (qi + q 2 ) 
It 

where 

q a = ^(1 + + v/i(l + £) 2 + {kl? + *)• (71) 

Therefore the prediction of the transient dispersion relation consists of two steps. 
Firstly the numerical resolution of Eq. d69t for a planar interface is performed by a 
Newton-Raphson method, yielding the transient front position This function 
is then introduced into Eq. dTOl) to obtain u(kl, t) in this adiabatic approximation. 

4.1.2 Results from the phase-field model 

We will use a PFM introduced by Losert et al. 1 109 1 for symmetric directional 
alloy solidification with constant miscibility gap. This model can be obtained from 
that of Sec. I2.3l bv making = (i.e. no kinetic dynamics) in Eq. |49] Once the 
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Figure 1 : Evolution of the interface position during the transient for (a) It /I = 12.5 
and (b) It 1 1 = 2.5, and convergence to the sharp-interface as e — > 0. 

thermal gradient is included, the equations for the phase field and the concentration 
take the following form in reduced variables: 



a i = an d a>2 = 0.6267 are integral constants obtained when performing 
the thin-interface limit [47 1. Note that the equation for the phase field evolution 
contains three parameters: two main control parameters coming from the sharp- 
interface model (It /I and do /I), and the model-specific parameter e. 

We will consider the transient from the rest, i.e. for an initial condition at to = 
consisting of an equilibrium solid-liquid planar interface located at £(i = 0) = 
l T II. In this case d{f,0) = Coo - AC (ui(r,0) = -1), and C 2 (r,0) = C M 
(u2(r, 0) = 0).That corresponds to taking u(f, 0) = — 1 as initial condition for 
the phase-field model. Trying to mimic a real experiment, we consider here that at 
t = the pulling velocity suddenly takes the final value V (1 in scaled variables). 

We perform numerical integration of the PFM Eqs. (1721731) with an explicit 
finite-differences scheme with Ax = 0.8 and At = 0.08 . We first study the 
ID dynamical evolution during the transient, comparing the results with the sharp- 
interface predictions. Fig. 1 presents the front position for two different values of 
the control parameter It /I = 12.5 and 2.5. For each case, we compare simulations 
for three different values of e (e = 0.5, e = 0.25 and e = 0.125) with the front 
position obtained with direct resolution of the integral equation (169b . Convergence 
to the sharp-interface limit can be observed as e decreases, and good agreement is 
found for a value of e = 0.125. 

We now estimate the transient dispersion relation u(kl,t) from PFM simula- 
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Figure 2: a) Localization of the linear growth regime of the mode kl = 1.25. b) 
Time evolution of two different modes kl = 1.25 and kl = 2.5 for It/I = 12.5. 

tions and compare it with the sharp interface prediction of Eqs. dTUl ITTl . To this 
end we simulate for each desired value of t a (planar) ID interface evolving from 
to = to t. At that moment we introduce a sinusoidal interface perturbation with 
wavevector kl, and continue the simulation in 2D. The spectral analysis of the front 
allows us to locate the regime where the mode evolution is linear. This is repre- 
sented in Fig. 2a, where three definite regions can be observed. From the linear 
region the value of the transient growth rate u(kl, t) is calculated. Fig. 2b shows 
the growth rate at different times obtained for two different modes kl = 1.25 and 
kl = 2.5 in the case of It/I = 12.5 and d^/l = 0.06923. Quantitative agreement 
is observed for all times. 

We can now introduce internal fluctuations in the diffusion equation of the 
phase-field by introducing the noises appearing in Eqs. (110113b . In simulations the 
growth of a range of wavelengths can be observed until a selected mode dominates. 
The spectral analysis of the front reveals quantitative agreement during the linear 
regime between the growth of each mode and the transient dispersion relation pre- 
dicted by the sharp interface model. In Fig. 3 we present the amplitude of one of 
the modes (kl = 2.5, It/I = 12.5) averaged for 15 different noise realizations. 
There is a good agreement in the regime with positive oj (i.e. were an increasing 
Ak(t) is predicted). 

4.2 Dendritic sidebranching in free solidification 

In the previous Section our interest was the initial destabilization of the planar 
interface in solidification experiments, and the influence of the internal noise. That 
transient has relevance in the selection of the wavelength of the dendritic array. 
In the present Section our interest will lie in some characteristics of the dendrite 
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Figure 3: Time evolution of the mode kl = 2.5 (It/I = 12.5) with the phase-field 
with noise, and comparison with the sharp-interface prediction. 



itself, in particular in those of the sidebranching. We study the noise-induced 
sidebranching scenario in the linear regime and in particular we focus on the issue 
of the internal vs. external origin of the noise [38]. 

In a frame of reference moving with the tip of the dendrite, sidebranching can 
be seen as a wave that propagates along the dendrite away from the tip at the tip's 
velocity. Therefore, an appropriate characterization is provided by the amplitude 
of the wave. Brener and Temkin [ 105 1 pointed out that experimentally observed 
sidebranching could be explained by considering noise of thermal origin. They 
analytically found that the growth of sidebranching amplitude in the linear regime 
behaves exponentially as a function of z 2 / 5 . Dougherty et al. II 101 studied side- 
branching in experiments with ammonium bromide dendrites. The sidebranching 
amplitude was found to increase exponentially up to a certain value of z, from 
which the linear theory is presumably no longer valid. It was also observed that 
side branches separated by more than about six times the mean wavelength are 
uncorrected. 

Bisang and Bilgram [ 1 1 1 1 found a quantitative agreement between the predic- 
tions for the linear regime in Ref. [ 105 1 and their results in experiments with xenon 
dendrites. They concluded that Brener and Temkin calculations describe correctly 
the sidebranching behaviour of dendrites for any pure substance with cubic sym- 
metry and thus thermal noise was the origin of the sidebranching observed in their 
experiments. 

Karma and Rappel [ 39 1 included thermal noise in a two-dimensional PFM for 
solidification controlled by heat diffusion, and obtained a good quantitative agree- 
ment between the computed sidebranching amplitude and wavelength as a function 
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of distance to the tip and the predictions of the linear WKB theory for anisotropic 
crystals in two dimensions. For a needle crystal shape x ~ z 3 / 5 , sidebranching 
amplitude was found to increase exponentially as a function of z 2 / 5 . 

However, the extension of the theory of dendritic sidebranching to the growth 
controlled by solutal diffusion shows that there are indications that in some exper- 
iments internal thermodynamical fluctuations could not account for the observed 
sidebranching activity [38 1. In that case, some other source of fluctuations, of ex- 
ternal origin, should be invoked. We study some of the consequences derived from 
adding a non-conserved noise source into a two-dimensional PFM. This noise is 
of very different nature than what one should employ to account for internal fluc- 
tuations. 

4.2.1 Sidebranching characteristics 

We have performed simulations of dendritic growth by employing a PFM for so- 
lidification with anisotropy included in the surface tension l20l . which has been 
taken as a = <x(0)(l + ^ su cos(49)). Whereas internal fluctuations would have ap- 
peared as a stochastic current (and, therefore, as a conserved term) in the equation 
for the diffusion field and an additional stochastic term in the phase field equation , 
like in Eqs. (I8l9t . we have chosen to model an external source of fluctuations by a 
non-conserved random term added to the equation for u. This random term reads 
simply / • r, where I denotes the amplitude of the noise, and r is a uniform random 
number in the interval [—0.5,0.5]. More details about the employed numerical 
procedure and the selection of parameters can be found in Ref. (38). 

A growing dendrite obtained by means of the PFM is shown in Fig. |4] Side 
branches appear at both sides of the main dendrite, yielding approximately a 90° 
angle with it, like it was observed in Ref. II 101 . Far down the tip one can clearly 
observe competition between branches which gives rise to a coarsening effect. In 
order to study the sidebranching induced by noise we have measured the half- 
width h z (t) of the dendrite at various distances z behind the tip as a function of 
time. The half-width of the dendrite as a function of time and its power spectrum 
P z (f) are shown in Fig. |5]for a distance z = 100 grid points. The data used to 
compute the power spectrums were six times the shown lengths. We also computed 
the cross-correlation function 

C{t') =< [h Lz (t + t') - h Lz ][h Rz (t) - h Rz ] > /a L a R , (74) 

where /iL 2 (i), h Rz (t), and a R are the half-width functions and their standard 
deviations for the two sides of the dendrite at the same distance from the tip. We 
found that C(0) is around 0.4 for points very close to the tip and that its value 
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Figure 4: Dendrite obtained with the PFM with an external noise term. 



decreases very quickly to when increasing z. The same behavior was observed 
in the experiments in Ref. II I OH . 

In Fig. |6]we show the logarithm of the square root of the area under the spectral 
peak as a function of z. This representation gives the behavior of the amplitude of 
the sidebranching as a function of the distance to the tip. The amplitude has a dif- 
ferent behavior for z up to 100 (linear regime) and for values greater than this one 
(nonlinear regime). In the region near the tip the amplitude increases exponentially 
with z 2 / 5 , while far down from the tip the growth rate of the amplitude is slower. 
Thus, the behavior of the data obtained in the simulations is consistent with the 
linear analysis carried out in Refs. l39l[T0 5 1 up to a certain value of z. 

In conclusion, despite thermal noise could not always be the main origin of 
sidebranching, its qualitative characteristics are common for noise-induced side- 
branching independently of its origin. 

4.3 Mesophase transitions in liquid crystals 

The growth of liquid crystal interfaces in transitions between mesophases is in 
many aspects analogous to that of solidification interfaces. The basic description is 
expected to be the same, lying the main differences in the parameter ranges, which 
often makes the liquid crystal case particularly suitable from an experimental point 
of view. From a theoretical point of view, a significant parameter difference with 
respect to the solidification case is the one associated to the diffusion coefficients, 
which are of the same order of magnitude in the two liquid crystal phases. An- 
other distinct feature is the presence of strong anisotropics, both at the interface, 
which can be faceted, and in the transport properties in the bulk. In this Section 
we will show how these anisotropics, introduced into a PFM, can account for the 
morphologies obtained in different growth conditions. 
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Figure 6: Square root of the area under the spectral peak as a function of ; 



28 







0.8 



Figure 7: CCH3. N-SmB interface (points with <fi = 0.5) obtained from simulation at 
subsequent times. A = 0.05. 



4.3.1 Anisotropic surface tension 

We have considered a liquid crystalline substance, CCH3, which presents a N- 
SmB phase transition at temperature T^s- The pattern formation during the N- 
SmB phase transition can be observed experimentally in a quasi two-dimensional 
geometry [29 30 1. The sample was initially set above the phase transition temper- 
ature T/vs, and was suddenly undercooled below this temperature in such a way 
that small smectic-B germs nucleated and started to grow. 

From the equilibrium shape one can obtain the corresponding angular depen- 
dence of the normalized surface tension r}{9) = (being 9 = the direction 
parallel to the smectic layers) by means of the Wulff-construction 11121 . The func- 
tion for CCH3 is HUE): r)(9) = 1.00 - 0.359 2 + 0.01(9 4 in the range |0| < vr/2, 
while 7](9) = rj(jr — 9). Clearly rj(9) has cusps at 9 = ±tt/2. 

Out of equilibrium, CCH3 presents three different morphologies, ranging from 
a slightly deformation of the equilibrium shape, to a butterfly-like morphology and, 
finally in the largest undercooling regime, to a four-fold dendritic shape ll29l l30l 

El. 

The nematic - smectic-B (N-SmB) phase transition is of first order, therefore 
standard models for solidification can be applied. Like in the previous section, we 
have employed the anisotropic PFM of Ref. [20|. In simulations the kinetic term 
was taken isotropic. The employed numerical procedure was described in [29|. 
We have simulated one quarter of the full experimental system by locating the 
initial smectic-B seed in the lower left corner. Symmetrical (reflecting) boundary 
conditions for <p and u have been imposed on the four sides of the system. 

The four main morphologies of CCH3 have been computationally reproduced 
on a qualitative level. At small dimensionless undercooling, simulations show a 
very slow growth of the germ, which at the observed times maintains a rectangle- 
like shape very similar to the equilibrium one in experiments, i.e., with two parallel 
facets and two rough convex sides (Fig. 0. The growth velocity of the facets is 
smaller than that of the convex parts of the interface. At slightly larger undercool- 
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Figure 8: CCH3. N-SmB interfaces ((f> = 0.5) obtained from simulation at subsequent 
times. A = 0.09. 
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Figure 9: CCH3. N-SmB interfaces ((f> = 0.5) obtained from simulation at subsequent 
times. A = 0.2. 



ings, the short sides undergo a first instability from convex to concave (Fig. At 
larger values of the undercooling the faceted sides start to bend adopting a slightly 
concave curvature and macroscopic facets disappear with the four corners opened 
up, forming a butterfly-like shape (Fig. |9). Finally, a completely developed four- 
fold dendrite is obtained in the large undercooling regime (Fig. ITOt . 

4.3.2 Anisotropic heat diffusion 

It has been proved in the previous Section the important role played by the surface 
tension anisotropy in the selection of the steady-state shape of a growing inter- 
face. Aside from this one and kinetic anisotropy, systems such as liquid crystals do 
present considerable anisotropics in their transport coefficients 1111311114111 15111161 
and their effects are potentially significant in interfacial growth phenomena. We 
address here the growth of dendrites subject to anisotropy in the heat diffusion 
coefficient, and its consequences on the growth morphologies [33]. 

When the anisotropic thermal diffusion is implemented in a PFM, the ^equation 
does not change and the u-equation should be slightly modified to consider a dif- 
fusion matrix Dij. The numerical details can be found in Refs. I33ll32l . 

The experimental system used to study the thermal diffusion anisotropy is the 
liquid crystalline substance CCH4, and in particular its N-SmB phase transition. 
At large undercoolings, some SmB seeds whose director directions do not match 
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that of the sorrounding N phase, present a non-reflection symmetry in their growth. 

With the modified PFM, it is possible to reproduce such experimental mor- 
phologies. This is observed in the numerical simulation (Fig. [TTl . Although 
the used set of parameters could differ from those real material ones, the quali- 
tative resemblance with experimental results is remarkably good [33 1. Axis x,y 
has been placed along the principal directions of the diffusion matrix, in order to 
have only diagonal terms in the equation. Anisotropy in the diffusion, defined 
as [Pmax — D m i n ) / D m i n , has taken a value of 0.5, very similar to experimen- 
tal estimates. The surface free energy has been rotated in the (x,y) plane, ac- 
counting for different orientations of the director in each phase. It can be derived 
from Fig. [^that the reflection symmetry has been broken by the inclusion of 
the anisotropic heat diffusion (previous simulations with rotated surface tension 
function and isotropic heat diffusion did not show the asymmetry in the growth 
velocities). In both experiments and simulations the most developed branches are 
systematically those growing in the direction of lowest heat diffusion. Therefore 
one can conclude that heat diffusion anisotropy favours dendritic growth in the 
lowest diffusion directions. This means that the relevant heat diffusion process is 
the one that occurs in the direction perpendicular to the axis of the dendrite. 



31 



Figure 1 1 : Simulation for the non-reflection symmetry. 



4.4 Results on Viscous Fingering 

The puipose of the present Section is to illustrate the practical use of the PFM in the 
case of viscous fingering, with direct quantitative tests of the thin-interface approx- 
imation to assess the usefulness of the approach in this case. Indeed, the fact that 
the model has the correct sharp-interface limit does not guarantee its practical use- 
fulness for several reasons. On the one hand, the stability of both the bulk phases 
and the kink profile must be assured, since this might not be the case in general. 
On the other hand, a direct empirical test is necessary to determine quantitatively 
how close a finite e situation is to the sharp-interface limit. This means finding a 
set of explicit quantitative criteria to choose all the nonphysical parameters in or- 
der to ensure a desired accuracy. Finally, it is interesting to find to what extent that 
model can provide quantitative results with reasonable computing efforts in actual 
simulations. 

4.4.1 Linear dispersion relation 

The first situation in which the model can be tested is the linear regime of a 
perturbed planar interface. The linear dispersion relation has been computed for 
vanishing viscosity contrast (c = 0). The sharp-interface model predicts a linear 
growth that does not depend on c. However, the phase field model should exhibit 
some dependence in the viscosity contrast related to the finite-e and -e corrections 
(see Sec. IV in Ref. B51 l. 

We use a single mode occupying the whole channel width (i.e., of wavelength 1 
and wavevector k = 2it) and then vary the dimensonless surface tension B in order 
to change the growth rate of that mode according to the Hele-Shaw dispersion 



32 




-1.0 



0.0 



0.5 



1.0 



Figure 12: Linear dispersion relation in scaled variables 



relation 



u = \k\(l-Bk 2 ). 



(75) 



This is physically completely equivalent to fixing the surface tension and varying 
the wavevector of the mode, as can be seen through the rescaling k' = VBk, uj' = 
\fBu>o = I k' | (1 — k' 2 ). However, it is numerically more convenient, since it allows 
one to use the same value of e for all the modes, because k is fixed. Thus, according 
to the finite-e and -e dispersion relation derived in Sec. IV of Ref. [43 1, 



We have used e = 0.01, e = 0.1 for the deviation from the sharp-interface one, ujq, 
to keep below a 10%. 

In Fig. El we present the linear dispersion relation thus obtained in scaled 
variables. The points (+) correspond to the measured growth rates for roughly a 
decade in the amplitude. Their deviation from the sharp-interface result of Eq. 
d75t (solid line) keeps below the desired 10% error and is fairly well quantitatively 
predicted by the thin-interface dispersion relation of Eq. dTBT i (dotted line). This 
quantitative agreement between theory and numerics is quite remarkable if we take 
into account that the thin-interface model is based on an asymptotic expansion in 
e. This good agreement is indeed an indication that the value of e used is in the 
asymptotic regime of the sharp-interface limit, as we will see more clearly in Fig. 



The growth rate values shown in Fig. El could still be refined by further de- 
creasing e and e. This is not only a theoretical possibility, but can also be done 
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e\k\V2l) + 0(c 2 ) + 0(e 2 ). 



(76) 
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Figure 13: Convergence to the Hele-Shaw value for the maximum of the linear dispersion 
relation curve 



in practice as we show in Fig. although the computation time increases as ex- 
plained above. Here, we study the convergence of the growth rate to (y-axis) for the 
maximum of uj'(k') (B = 8.443... x 1(T 3 ) to the Hele-Shaw result to = 4.188... 
(left upper corner) as we decrease e (x-axis) and e (various symbols). The empty 
symbols have been obtained with Ax = e, whereas the filled ones correspond to 
Ax = e/2. The growth rates obtained with Ax = e are always below the ones 
for Ax = e/2, probably because of the stabilizing effect of the mesh size. If Ax 
is further decreased the differences with the values computed with Ax = e/2 are 
tiny, whereas the gap between the Ax = e and the Ax = e/2 points is pretty large 
(clearly more than the differences between distinct symbols — distinct e values — 
or adjacent values of e). This means that the discretization has practically con- 
verged to the continuum model for Ax = e/2, but not for Ax = e. That is the 
reason why we have used Ax = e/2 in Fig. ^] 

Moreover, the Ax = e/2 points should be described by the thin-interface 
model of Eq. (f76l . and this is indeed the case for small enough values of e. To 
visualize this we have plotted the thin-interface prediction for e = (dashed line), 
which is, of course, a straight line in e. Each set of points with a same e value 
clearly tends to align parallel to this line as e decreases — as Eq. dToT i predicts — , 
whereas they curve up and even cross the line for large values of e, for which we 
are beyond the asymptotic regime of validity of Eq. (f76t . apparently ending near 
e = 0.01. This makes e = 0.01 very suitable for simulations, and confirms it to 
be within the asymptotic regime as we pointed out above. Note as well how the 
growth rate increases with decreasing values of e within the same value of e (ver- 
tical columns of points), and how it approaches the dashed line in good agreement 
with the values predicted by Eq. dToT i for values of e within the asymptotic regime. 
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Figure 14: Time evolution of the interface a) c = b) c = 



.8 



4.4.2 Finger dynamics 

The quantitative test in the nonlinear regime are not so simple due to the lack of 
exact analytical solutions. For the case of asymptotic Saff man-Taylor fingers, for 
instance, only partial numerical information is available, in particular for varying 
viscosity contrast. A detailed comparison with existing evidence was carried out 
in Ref . |H 

At a qualitative level, however, it is illustrative to plot an example of multifinger 
dynamics, in particular to emphasize the importance to have the viscosity contrast 
as an explicit parameter in the model. Indeed, the scenario of finger competition in 
the high viscosity contrast limit, where larger fingers screen out the smaller ones, 
is actually valid only for c very close to 1. In contrast, the persistence of the small 
fingers, which is characteristic of low viscosity contrasts, seems to be the generic 
case fTOfl . 
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Here we use a somewhat experimentally realistic initial condition consisting of 
a superposition of sinusoidal modes with random, uniformly distributed amplitudes 
between -0.005 and +0.005 for each wavelength A = 1, |, — i.e., in the linear 
regime and random phases. For B = 10~ 3 the most amplified of these wavelengths 
will be A = jj, so that we expect 3 unequal fingers to appear and there is a chance 
for mode interaction and competition to set in. Wavelengths below A = 0.161... 
are stable and will decay. We include some of them anyway. Then all modes are 
added up to find the interface position. The stream function predicted by the linear 
theory is also obtained by adding up the stream function of each mode, but all with 
their peaks centered at the same final interface position, to avoid the formation of 
more than one peak of the stream function across the interface. 

Since harmonics of the channel width are present, we have to refine the e used 
in the linear regime. We use e = 0.00625, with Ax = e to save computation time. 
The value of e is quite crude (e = 0.5) for the equal viscosities case c = 0, Fig. 
[14k . and especially for the high viscosity contrast (c = 0.8) run Fig. IT4b (e = 0.2). 

The results are shown in Figs. [T4h.b for 0.15 < t < 1.25 at constant time 
intervals 0.1 (dots). The last interface is emphasized in bigger points (+) and the 
solid line corresponds to the initial condition for t = 0. As we can see, the initial 
condition happens to have six maxima. Rather quickly, only three of them are left 
out as predicted by linear stability, even before entering the non-linear regime, in 
which these maxima elongate into well developed fingers. For vanishing viscosity 
contrast (c = 0, Fig. UAk). there is no apparent competition, in agreement with 
experimental [ 100 1 and numerical l90ll95lfT0Tl evidence. Longer and shorter fin- 
gers all advance. Shorter fingers might not advance so quickly, but they expand to 
the sides, so they clearly keep growing. In contrast, starting with the same initial 
condition that for c = 0, the c = 0.8 run (Fig. Et>) shows competition between 
fingers of the less viscous fluid advancing into the more viscous one, as it is known 
to happen in the Saff man-Taylor problem. The shorter finger now also expands 
laterally, but it soon begins to move backwards as a whole. Longer times may lead 
to the pinch-off of droplets in both cases 11011 . 

In summary the basic criteria to control the closeness to the sharp interface 
limit are en « 1, nsp? << 1- More precisely, we find numerically that the 
thin-interface model is accurate in the linear regime with an error below below 
10% if one satisfies the conditions ek < 0.06, nj^p — 0-016. The method 
could be made more efficient by using an adaptive mesh or (possibly) by can- 
celling out the corrections to the sharp-interface equations remaining (i.e., other 
than the Allen-Cahn law) in the thin-interface model. For high viscosity contrasts, 
c ~ 1, a distinct model could possibly be more efficient. From the recent progress 
in boundary-integral methods [98], however, the use of PFM may not be signifi- 
cantly advantageous in this problem except maybe for the natural incorporation of 



36 



interface pinch-off. Nevertheless, the application of the PFM approach in three- 
dimensional viscous flows in porous media, is yet an unexplored and promising 
future direction. 

5 Conclusions 

Phase-field modeling of interface dynamics in nonequilibrium systems has proved 
to be a very useful tool. In this chapter we have reviewed different aspects of 
this active field of research. Phase-field models introduce an auxiliary field to al- 
low for a diffuse interface. This additional field is coupled to the physical fields 
in such a way that the macroscopic, sharp-interface boundary conditions are ef- 
fectively reproduced. In this way, a free-boundary problem is transformed into a 
set of partial differential equations, usually much simpler to handle. In this chap- 
ter a phase-field model for solidification has been explicitly derived from a free 
energy. As a second example, a phase-field model for the Saff man-Taylor prob- 
lem has been presented. In this case the derivation is not based on a free-energy 
functional. The thin-interface approximation, which includes corrections on the 
interface thickness, have been discussed for these two examples. We have also pre- 
sented some illustrative applications such as the study of fluctuations in directional 
solidification of alloys and dendritic sidebranching in free solidification, growth of 
mesophases in liquid crystals and the dynamics of viscous fingering. 

At the present level of development of phase-field model techniques, these have 
been shown to be quantitatively accurate with reasonable computational cost, and 
usually advantageous with respect to other techniques in different aspects. Phase- 
field models are simple to be implemented in a computer and contain a complete 
description of the full nonlinear and nonlocal properties of the macroscopic free- 
boundary problems. In addition, in the phase-field framework it is relatively simple 
to introduce additional perturbations, such as fluctuations and static disorder, or 
modifications of the boundary conditions. Future lines of research will certainly 
include quantitative studies of interfaces in three dimensions, the consideration of 
complex fluids and biological systems, and possibly applications at the nanoscale. 
Undoubtly the phase-field approach holds the promise of fruitful applications to 
new problems in years to come. 
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